/*
 *  Vector.cpp
 *  T3nsors
 *
 *  Created by Michael Barriault on 09-10-07.
 *  Copyright 2010 University of Guelph. All rights reserved.
 *
 */

#include "Vector.h"
#include "Macros.h"
#include "Others.h"
#include <cstdarg>
#include <cstdio>
#include <cmath>
using std::string;

Vector::Vector(void) {
	return;
}

Vector::Vector(Domain iO, int iN, std::string iid) {
	O = iO;
	N = iN;
	id = iid;
	T.assign(N, Scalar(O,id));
	FOR(n,N)
		T[n].id = T[n].id + itos(n);
}

Vector::~Vector(void) {
	T.clear();
}

Scalar& Vector::operator()(int i) {
	return T[i];
}

double Vector::rad(int o) {
	double r = 0;
	FOR(n,N) r += T[n][o]*T[n][o];
	return sqrt(r);
}

double& Vector::operator()(int m, int a, ...) {
	va_list va;
	va_start(va, a);
	int t=0;
	FOR(i,O.N) {
		int A = a;
		for ( int j=i-1; j>=0; j-- )
			A *= O[j];
		t += A;
		a = va_arg(va, int);
	}
	va_end(va);
	return T[m][t];
}

Vector operator-(Vector x) {
	Vector y(x.O,x.N,(string)"-"+x.id);
	FOR(i,x.N)
		y(i) = -x(i);
	return y;
}

Vector operator+(Vector x, Vector y) {
	CheckDims(x.O,y.O); CheckDims(x.N, y.N);
	Vector z(x.O,x.N,x.id+"+"+y.id);
	FOR(i,x.N)
		z(i) = x(i) + y(i);
	return z;
}

Vector operator-(Vector x, Vector y) {
	CheckDims(x.O,y.O); CheckDims(x.N, y.N);
	Vector z(x.O,x.N,x.id+"-"+y.id);
	FOR(i,x.N)
		z(i) = x(i) - y(i);
	return z;
}

Scalar operator*(Vector x, Vector y) {
	CheckDims(x.O, y.O); CheckDims(x.N, y.N);
	Scalar z(x.O);
	FOR(i,x.N)
		z = z + x(i)*y(i);
	z.id = x.id+"*"+y.id;
	return z;
}

Vector operator%(Vector x, Vector y) {
	CheckDims(x.O, y.O); CheckDims(x.N, y.N);
	Vector z(x.O, x.N, x.id+"%"+y.id);
	z(0) = x(1)*y(2)-x(2)*y(1);
	z(1) = x(2)*y(0)-x(0)*y(2);
	z(2) = x(0)*y(1)-x(1)*y(0);
	return z;
}

Vector operator*(Vector x, Scalar y) {
	CheckDims(x.O, y.O);
	Vector z(x.O,x.N);
	FOR(i,x.N)
		z(i) = x(i)*y;
	return z;
}

Vector operator*(Scalar x, Vector y) {
	return y*x;
}

Vector operator*(Vector x, double y) {
	Vector z(x.O,x.N);
	FOR(i,x.N)
		z(i) = x(i)*y;
	return z;
}

Vector operator*(double x, Vector y) {
	return y*x;
}

Vector operator/(Vector x, Scalar y) {
	CheckDims(x.O, y.O);
	Vector z(x.O,x.N);
	FOR(i,x.N)
		z(i) = x(i)/y;
	return z;
}

Vector operator/(Vector x, double y) {
	return x*(1./y);
}

Vector diffOld(Scalar x) {
	Vector dx(x.O, x.O.N, (std::string)"d"+x.id);
	FOR(p,dx.N) {
		int dg = x.O(p);
		int dc = x.O(p-1);
		int G = x.O()/dg;
		int C = x.O[p];
		FOR(g,G) {
			FOR(m,dc) {
				int i = g*dg + m;
				if ( x.O.p(p) )
					dx(p)[i] = (x[i+dc]-x[i+(C-1)*dc])/(2*dx.O.d(p));
				else
					dx(p)[i] = (x[i+dc]-x[i])/dx.O.d(p);
				for ( int c=1; c<C-1; c++ ) {
					i = g*dg + m + c*dc;
					dx(p)[i] = (x[i+dc]-x[i-dc])/(2*dx.O.d(p));
				}
				i = g*dg + m + (C-1)*dc;
				if ( x.O.p(p) )
					dx(p)[i] = (x[i-(C-1)*dc]-x[i-dc])/(2*dx.O.d(p));
				else
					dx(p)[i] = (x[i]-x[i-dc])/dx.O.d(p);
			}
		}
	}
	return dx;
}

#define D1 0.6666666666666666666666
#define D2 -0.083333333333333333333
#define D00 -1.4117647058823529411764705882353
#define D01 1.7352941176470588235294117647059
#define D02 -0.23529411764705882352941176470588
#define D03 -0.088235294117647058823529411764706
#define D10 -0.5
#define D12 0.5
#define D20 0.093023255813953488372093023255814
#define D21 -0.68604651162790697674418604651163
#define D23 0.68604651162790697674418604651163
#define D24 -0.093023255813953488372093023255814
#define D30 0.030612244897959183673469387755102
#define D32 -0.60204081632653061224489795918367
#define D34 0.65306122448979591836734693877551
#define D35 -0.081632653061224489795918367346939

Vector diffNew(Scalar x) {
	Vector dx(x.O, x.O.N, (std::string)"d"+x.id);
	FOR(p,dx.N) {
		int dg = x.O(p);
		int dc = x.O(p-1);
		int G = x.O()/dg;
		int C = x.O[p];
		FOR(g,G) {
			FOR(m,dc) {
				int i = g*dg + m;
				dx(p)[i] = (D00*x[i]+D01*x[i+dc]+D02*x[i+2*dc]+D03*x[i+3*dc])/dx.O.d(p);
				dx(p)[i+dc] = (D10*x[i]+D12*x[i+2*dc])/dx.O.d(p);
				dx(p)[i+2*dc] = (D20*x[i]+D21*x[i+dc]+D23*x[i+3*dc]+D24*x[i+4*dc])/dx.O.d(p);
				dx(p)[i+3*dc] = (D30*x[i]+D32*x[i+2*dc]+D34*x[i+4*dc]+D35*x[i+5*dc])/dx.O.d(p);
				for ( int c=4; c<C-4; c++ ) {
					i = g*dg + m + c*dc;
					dx(p)[i] = (-D2*x[i-2*dc]-D1*x[i-dc]+D1*x[i+dc]+D2*x[i+2*dc])/dx.O.d(p);
				}
				i = g*dg + m + (C-1)*dc;
				dx(p)[i-3*dc] = -(D35*x[i-5*dc]+D34*x[i-4*dc]+D32*x[i-2*dc]+D30*x[i])/dx.O.d(p);
				dx(p)[i-2*dc] = -(D24*x[i-4*dc]+D23*x[i-3*dc]+D21*x[i-dc]+D20*x[i])/dx.O.d(p);
				dx(p)[i-dc] = -(D12*x[i-2*dc]+D10*x[i])/dx.O.d(p);
				dx(p)[i] = -(D03*x[i-3*dc]+D02*x[i-2*dc]+D01*x[i-dc]+D00*x[i])/dx.O.d(p);
			}
		}
	}
	return dx;
}

Vector diff(Scalar x) {
	bool runNew = true;
	FOR(n,x.O.N) if ( x.O[n] < 8 ) { runNew = false; break; };
	if ( runNew ) return diffNew(x);
	else return diffOld(x);
	return diffOld(x);
}

Scalar diss(Scalar x) {
	Vector dx(x.O, x.O.N, (std::string)"d4"+x.id);
	FOR(p,dx.N) {
		int dg = x.O(p);
		int dc = x.O(p-1);
		int G = x.O()/dg;
		int C = x.O[p];
		FOR(g,G) {
			FOR(m,dc) {
				int i = g*dg + m;
				dx(p)[i] = (2*x[i]-4*x[i+dc]+2*x[i+2*dc])/dx.O.d(p);
				dx(p)[i+dc] = (-2*x[i]+5*x[i+dc]-4*x[i+2*dc]+x[i+3*dc])/dx.O.d(p);
				for ( int c=2; c<C-2; c++ ) {
					i = g*dg + m + c*dc;
					dx(p)[i] = (x[i-2*dc]-4*x[i-dc]+6*x[i]-4*x[i+dc]+x[i+2*dc])/dx.O.d(p);
				}
				i = g*dg + m + (C-1)*dc;
				dx(p)[i-dc] = (x[i-3*dc]-4*x[i-2*dc]+5*x[i-dc]-2*x[i])/dx.O.d(p);
				dx(p)[i] = (2*x[i-2*dc]-4*x[i-dc]+2*x[i])/dx.O.d(p);
			}
		}
	}
	Scalar Dx(x.O, dx.id);
	FOR(p,dx.N) Dx = Dx + dx(p);
	return Dx;
}

Vector diss(Vector x) {
	Vector dx(x.O, x.O.N, (std::string)"d4"+x.id);
	FOR(p,dx.N) dx(p) = diss(x(p));
	return dx;
}

